// STL
#include <algorithm>               // std::min_element
#include <cmath>                   // std::sqrt, std::pow
#include <iterator>                // std::distance, std::begin, std::end
#include <utility>                 // std::pair, std::make_pair
#include <vector>                  // std::vector

// jScience
#include "jsvd.h"                  // svdfit()


static void func(double x, double afunc[], int ma);


//
// DESC: Fits the objective function to a function defined below, and returns the function and the minimum point on this function.
//
// INPUT:
//     std::vector<double> E       : Objective function
//     std::vector<double> stddevE : Standard deviation of the objective function
//
// OUTPUT:
//     std::vector<double>::size_type indx : Index of the minimum of f(x)
//
// NOTES:
//     ! Both this routine and func() below are *very* sloppy in the sense that the use C-style offsets, static arrays, etc.
//
inline std::pair<
                 std::vector<double>,
                 std::vector<double>::size_type
                 > fit_validation_error(
                                        const std::vector<double> E,
                                        const std::vector<double> stddevE
                                        )
{
    double x[E.size()+1];
    double y[E.size()+1];
    double sig[E.size()+1];

    int ma = 5;
    double a[ma+1];
    double chisq;

    double **u, **v;
    u = new double * [E.size()+1];
    v = new double * [ma+1];

    for( decltype(E.size()) i = 0; i < (E.size()+1); ++i )
    {
        u[i] = new double [ma+1];
    }

    for( int i = 0; i < (ma+1); ++i )
    {
        v[i] = new double [ma+1];
    }

//    double u[E.size()+1][ma+1];
//    double v[ma+1][ma+1];
    double w[ma+1];

    for( decltype(E.size()) i = 0; i < E.size(); ++i )
    {
        // starting x at 0 or 1 is irrelevant, as long as it is consistent with the evaluation of func() below
        x[i+1]   = static_cast<double>(i+1);
        y[i+1]   = E[i];
        sig[i+1] = stddevE[i];
    }

    svdfit( x, y, sig, E.size(), a, ma, u, v, w, &chisq, func );

    // << JMM: there is probably an analytical solution ... >>
    std::vector<double> f;
    for( decltype(E.size()) i = 0; i < E.size(); ++i )
    {
        double afunc[ma+1];
        func( static_cast<double>(i+1), afunc, ma );

        double fx = 0.0;
        for( int j = 0; j < ma; ++j )
        {
            fx += (a[j+1]*afunc[j+1]);
        }

        f.push_back(fx);
    }

    for( decltype(E.size()) i = 0; i < (E.size()+1); ++i )
    {
        delete [] u[i];
    }

    delete [] u;

    for( int i = 0; i < (ma+1); ++i )
    {
        delete [] v[i];
    }

    delete [] v;

    std::vector<double>::size_type indx = static_cast<std::vector<double>::size_type>( std::distance( std::begin(f), std::min_element( std::begin(f), std::end(f) ) ) );

    return std::make_pair(f, indx);
}

// f(x) = a + b*x^0.5 + c*x + d*x^1.5 + e*x^2
// << JMM: I don't like using so many parameters! >>
static void func(double x, double afunc[], int ma)
{
    afunc[1] = 1.0;
    afunc[2] = std::sqrt(x);
    afunc[3] = x;
    afunc[4] = std::pow(x, 1.5);
    afunc[5] = x*x;

//    afunc[1] = 1.0;
//    afunc[2] = x;
//    afunc[3] = x*x;
//    afunc[4] = x*x*x;
}


